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Abstract 

This paper presents a survey of ocean simulation and rendering methods in computer graphics. To model and 
animate the ocean 's surface, these methods mainly rely on two main approaches: on the one hand, those which 
approximate ocean dynamics with parametric, spectral or hybrid models and use empirical laws from oceano- 
graphic research. We will see that this type of methods essentially allows the simulation of ocean scenes in the 
deep water domain, without breaking waves. On the other hand, physically-based methods use Navier-Stokes 
Equations (NSE) to represent breaking waves and more generally ocean surface near the shore. We also describe 
ocean rendering methods in computer graphics, with a special interest in the simulation of phenomena such as 
foam and spray, and light's interaction with the ocean surface. 



Categories and Subject Descriptors (according to ACM 
CCS): 1.3.7 [Computer Graphics]: Three-Dimensional 
Graphics and Reahsm — Animation 1.3.8 [Computer Graph- 
ics]: Apphcations — 



1. Introduction 

The main goal of computer graphics is to reproduce in the 
most true-to-life possible way the perceived reality with all 
the complexity of natural phenomena surrounding us. The 
ocean's complexity is mainly due to a highly dynamic be- 
havior. Ranging from a quiet sea to an agitated ocean, from 
small turbulent waves to enormous shorebreaks, the dynamic 
motion of the ocean is influenced by multiple phenomena oc- 
curring at small and large scales. For several centuries, scien- 
tists have tried to understand and explain these mechanisms. 
In oceanographic research, physicists define the behavior of 
the ocean surface depending on its location: in deep ocean 
water areas (far from the coast), intermediate areas or shal- 
low water areas (close to the shore). This classification char- 
acterizes wave motions with different parameters. In deep 
waters, the free surface defined by the interface between air 
and water is generally subjected to a large oscillatory behav- 
ior, whereas in shallow waters waves break near the shore. 
Representing the visual complexity of these phenomena is a 



challenge, and the last 30 years have seen computer graphics 
evolve in order to address this issue. 

Different models can be used to represent ocean dynam- 
ics: parametric description, spectral description as well as 
models from Computational Fluid Dynamics (CFD) and 
more specifically Navier-Stokes Equations (NSE). The first 
category aims at computing the path of water particles and 
describes the free surface with parametric equations based 
on real observations, obtained from buoys or satellite mea- 
surements [Bie52]. The second category approximates the 
state of the sea by using waves spectrum [PM64, HBB*73, 
BGRV85] and computes waves distribution according to 
their amplitudes and frequencies. Finally, NSE can repre- 
sent dynamics of all types of fluid, including the dynamical 
behavior of the ocean. 

Ocean simulation methods in the computer graphics do- 
main can therefore be classified into two main categories: 
parametric/spectral methods that use oceanographic mod- 
els, and physically-based methods relying on NSE. Paramet- 
ric/spectral models work best in deep waters where they ac- 
curately represent the periodical motion of the sea. But since 
they do not take into account the interactions with the bottom 
of the sea in shallow waters, only physically-based methods 
deriving approximate solutions from NSE can reproduce the 
complexity of ocean dynamics near the shore. 
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Another important characteristics of oceans for computer 
graphics is their complex optical properties. For example, 
the color of ocean waters, varying from green to deep blue, is 
related to the concentration of phytoplankton particles. Sev- 
eral other phenomena, such as foam, sprays, water properties 
(turbidity, bubbles, . . . ) and light-water exchanges must also 
be simulated at the rendering stage. Addressing the simula- 
tion of these phenomena is a precondition to obtain realistic 
ocean scenes. 

This paper presents a survey of research works in ocean 
simulation and rendering in computer graphics. This in- 
cludes different methods specifically designed for ocean 
scenes, but also more general water simulation techniques 
that can be applied to ocean simulation. Papers presenting 
specific methods for fluid simulation, intended for exam- 
ple for rivers [TGOl, YNBH09] or fountains [WCHJ06] are 
not covered here; interested readers can also refer to [Igl04] 
where a more general survey of water simulation techniques 
is presented. 

In sections 2 and 3, we will focus on the methods specif- 
ically intended to model and animate the detailed surface 
of the ocean. We will use the classification presented in the 
previous paragraph: parametric/spectral methods describing 
the surface using models from oceanography, and NSE- 
based methods that simulate the dynamic behaviour of ocean 
waves. Section 4 is dedicated to realistic ocean rendering, 
particularly the representation of foam and sprays and the 
simulation of different light-water interactions. Finally, we 
will conclude by presenting possible perspectives for these 
techniques. 



2. Ocean dynamics simulation in deep water 

The methods presented in this section use theoretical models 
and/or experimental observations to describe the ocean sur- 
face in deep water and represent swell effects. We can sub- 
divide these methods into three main categories: first, those 
describing the ocean surface directly in the spatial domain, 
then those describing the surface in spectral domain and fi- 
nally hybrid methods combining the two. As we will show in 
the following section, spatial methods use a height map com- 
puted as a sum of periodical functions, and animated with a 
simple phase difference, in order to represent the ocean sur- 
face. Spectral approaches use a wave spectrum to describe 
the surface in the spectral domain and a Fourier Transform 
to obtain its transformation in the spatial domain. Finally 
the combination of the two, called hybrid methods, produces 
convincing surfaces that can be animated easily. 



2.1. Spatial domain approaches 

The main goal of spatial domain approaches is to represent 
the geometry of the water surface using a sum of periodical 
functions evolving temporally using a phase difference. 




Figure 1: Ocean surface obtained in [FR86] 




Figure 2: Shapes of waves obtained using Eq. 1 



2.1.1. Early works 

This idea to combine series of sinusoids with high and low 
amplitudes was first proposed by Max [Max81]. Ocean sur- 
face is represented as a height map with height y = h{x,z,t) 
computed at each point (x, z) at time t by: 



h{x, z,t) = -yo + Yj ^icosilci^x + ki_z - w,-f ) 
/■=! 



(1) 



where Nw is the total number of waves. A, is the amplitude of 
the ('-th wave, k, = {ki^_ , ki^ ) its wave vector, w; its pulsation 
and yo is the height of the free surface. The x-axis is oriented 
horizontally and points towards the coastline, the y-axis is 
vertical, and the z-axis is horizontal and aligned with the 
coastline. For each wave, the shape of the curve defined by 
the motion of a single point depends directly on the product 
between the amplitude A/ and the wave number A:,- = 1 1 fc/ 1 1 . If 
kiAi < 0.5, this path is similar to a trochoid. If A:,A, = 0.5, the 
shape is a cycloid. In all other cases (A:,A,- > 0.5), this path 
cannot represent a realistic motion (see Figure 2). 

In order to obtain a realistic effect, the wave vector fc,- of 
each wave is computed using scattering relationship in deep 
water domain, i.e. supposing that the bottom of the sea is 
located at infinite depth: 



(2) 
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with g the gravitation constant and L; the wavelength of each 
individual wave. 

The same idea was developed with a simple bump map- 
ping approach, where normal vectors on a planar mesh are 
transformed using a sum of 20 cycloids [SchSO]. However, 
assuming that the bottom of the sea is at infinite depth limits 
the use of these methods since they don't include more com- 
plex phenomena such as breaking waves or waves refraction 
near the shore. Peachey [Pea86] introduces a depth parame- 
ter to compute the wave vector ki of each wave, using Airy 
wave theory: 

ki= 271/ J ^tanh— (3) 
V 271 Lj 

where d is the depth of a point related to the bottom of the 
sea. 

Fournier and Reeves [FR86] suggested modifying Gerst- 
ner's theory of waves, where water particles positions are 
given by: 

I x = xo-Ae''^'<'sm{txQ-wt) 
I y = yo~ cos {kzo — wt) 

with X (respectively y) the horizontal (resp. vertical) coordi- 
nate of a water particle at time t, xq and yo its coordinates 
at rest, A the wave amplitude, k the wave number and w the 
wave pulsation. Fournier and Reeves enhance this model by 
taking into account the transformation of the path of water 
particles following the topological changes of the sea bed, 
and by transforming their circular path into a more realistic 
elliptic motion. This method permits to control the waves' 
shape, more or less crested, through the use of different pa- 
rameters, and therefore yields a more realistic result (see 
Figure 1). 

Gonzato and Le Saec [GS99] modify this model when the 
starting point of a plunging breaking wave is detected near 
the shoreline (i.e. if the wave's crest starts to curl over). This 
phenomenological modification results in the addition of two 
local functions applied to the wave's shape: a stretch func- 
tion imitates Biesel law by progressively stretching the wave 
along its crest, and a plunging function simulates gravity. 

Ts'o and Barsky [TB87] suggested that the refraction of 
waves could be represented by the principle of light refrac- 
tion formulated by Descartes. Their approach, called wave 
tracing, consists in generating a spline surface by casting 
rays from the skyline in a uniform 2D grid, progressing with 
Bresenham's algorithm. The deviation of a ray is computed 
according to the depth difference of a cell to the next. How- 
ever, a major drawback with this method is the lack of de- 
tails on the generated surface when rays diverge significantly 
from a straight line, since in that case the number of rays 
defining the surface is too low. Gonzato and Le Saec [GSOO] 
address this problem by generating new rays in undersam- 
pled areas, offering a better representation of the ocean sur- 
face around bays or islands for instance. This method called 



"Dynamic Wave Tracing" can handle wave reflexion and 
diffraction due to obstacles. This approach was also used 
by Gamito and Musgrave [GM02] to extract a phase map 
used in a parametric model, which can be animated eas- 
ily. The height map representing the ocean surface can be 
rendered using enhanced ray tracing [GSOO] or sphere trac- 
ing [DCG07] methods. 

2.1.2. GPU implementations 

Series of periodic functions are particularly well suited to 
GPU computations, and in the last few years research works 
describing real time implementations have emerged. 

Chen et al [CLW07] use four sinusoids, computed in a 
vertex shader and transferred to a bump map in a pixel shader 
to simulate ripples. Salgado and Conci [SCO?] describe a 
real-time implementation of Foumier's method [FR86] us- 
ing a vertex shader. Schneider and Westermann [SWOl] 
compute a displacement map using Perlin noise in a ver- 
tex shader on the GPU to obtain interactive results. Isidoro 
et al [IVB02] use a precomputed mesh perturbed by four 
sinusoids with low frequencies in a vertex shader. Ripples 
are obtained by using a bump map combining multiple tex- 
tures. When using a limited number of input functions, visu- 
ally convincing surfaces are obtained in interactive time (see 
Figure 3). Cui et al [CYCXW04] follow the same approach 
combined with an adaptive tessellation of the surface using 
the method of Hinsinger et al [HNC02] (described in Section 
2.3). 

Chou and Fu [CF07] described a new GPU implementa- 
tion for ocean simulation in order to take into account vari- 
ous interactions between the ocean surface and fixed or dy- 
namic objects. The aim of this approach is to combine a par- 
ticle system with a height-field obtained with a sum of si- 
nusoids computed at each vertex of the surface mesh. Parti- 
cles move according to the motion of each vertex but also to 
hydrodynamical forces [Bro9 1 , FLS04] due to solid objects 
interacting in the scene. 

2.1.3. Adaptive schemes 

In order to limit computation time, adaptive schemes have 
been proposed to limit computations only in visible areas 
and/or to reduce the number of periodic functions used to 
represent the surface dynamically, depending on the distance 
to the viewer. This approach consists in discarding high fre- 
quencies in distant areas, also limiting aliasing effects such 
as moire patterns visible on the horizon. The "grid project" 
concept proposed by Johanson [Joh04] consists in project- 
ing the pyramid of vision onto the height map representing 
the water surface. Distant areas are automatically sampled 
at low resolution, whereas high resolution is used near the 
viewer. 

Finch [Fin04] and Lee etal [LGL06a,LGL06b] propose a 
Level of Details (LOD) tessellation by adapting the resolu- 
tion of the surface to the distance between the viewer and the 
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Figure 3: Real-time ocean simulation from [KryOS] 



horizontal plane; using the angle between the camera and the 
surface, vertices outside the pyramid of vision are also dis- 
carded before rendering stage. The implementation by Lee et 
al provides a 40% gain in computation time compared to a 
fixed-resolution approach, and outperforms the implementa- 
tion of Johanson [Joh04] (tests were conducted on a Pentium 
4 at 3.0GHz and a GeForce 6600 GPU). 

Rryachko [KryOS] use vertex textures, by computing two 
height maps representing the global motion of waves and 
fine scale details. The resolution is based on the distance to 
the viewer, allowing for fine-tune computations (see Figure 
3). 



2.2. Fourier domain approaches 

This type of representation consists in describing the ocean 
surface using a spectral distribution of waves, obtained from 
theoretical or measured data. The spatial representation is 
then computed using the most significant frequency com- 
ponents and an inverse fast Fourier transform (IFFT), giv- 
ing a geometric characterization of the surface as a dynamic 
height map defined for each point X by: 



h{X,t) = '£h{k,t)e' 
1=1 



iliX 



(5) 



where A^, is the number of spectral components, k is the wave 
vector and li{k,t) is the amplitude of the Fourier component 
obtained from a theoretic wave spectrum. 

2.2.1. General metliods 

This idea was proposed by Mastin et al [MWM87] with 
the Pierson-Moskowitz spectrum [PM64]. The amplitude of 
each Fourier component at time / = is given by: 



hoik) ■ 



167l4/5 



(6) 



where a = 0.0081 is called Phillips' constant, g is the grav 
itation constant. The frequency peak /„, is defined by: 

0.1 3g 



fm — ■ 



(7) 



with Uio the wind speed measured at 10 meters above the 
surface. 

Premoze and Ashikhmin [PAOO] follow the same idea 
by using the JONSWAP spectrum [HBB*73], a modified 
version of Pierson-Moskowitz's which reduces frequencies 
and thus raises waves' amplitudes. The amplitude of each 
Fourier component at time t = in that case is slightly dif- 
ferent: 



-if-fm 



hoik) ■ 



167l4/5 



ln(y)e 



(8) 



with: 



0.07 if/ 

0.09if/>/m 



The method presented by Tessendorf [TesOl] was notably 
used in famous computer-generated scenes for Waterworld 
and Titanic productions. Again, this approach uses spectral 
components to describe the surface, computed by a Gaussian 
pseudo-random generator and a theoretic wave spectrum due 
to Phillips: 



Phik) 



(9) 



with Yr and Y; defined as constants. The spectrum Ph{k) is 
given by: 



Ph{k) = C 



(10) 



where C is a constant, k is the wave vector, k is the wave 

V- 

number, w is the wind direction, V is its speed and L = — . 
The choice of a random generator is justified by the Gaussian 
distribution of waves often observed in deep sea areas. In 
order to obtain fine details on the surface, each component 
is slightly perturbed by a noise function, producing more or 
less crested waves at any resolution (see Figure 4). Cieutat et 
al [CGG03] enhance this model in order to take in account 
water waves forces applied to a ship, implemented in a ship 
training simulator. 

Although results obtained with this approach seem visu- 
ally convincing for a distant viewer, when zooming towards 
the surface, waves' shapes tend to look unrealistic since the 
resolution of the heightmap is fixed. In that case it becomes 
necessary to generate a new surface using more spectral 
components. This problem is addressed by interactive meth- 
ods proposed recently. 

2.2.2. Level-Of -Detail and GPU implementations 
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Figure 4: Ocean surface obtained in [TesOl ] 



Hu et al [HVT*06] apply a LOD approach to Tessendorf's 
method by representing the ocean with two surfaces. The 
first one, with a high, fixed resolution, is animated using a 
displacement map stored in a vertex shader, while the second 
one is sampled adaptively and applied as a bump map stored 
in a pixel shader. This approach produces a large ocean sur- 
face in real-time by animating only visible parts (at approx. 
100 frames/sec. on GeForce 3 GPU card). Mitchell [Mit05] 
implemented Tessendorf's method on GPU by only consid- 
ering frequencies generating a significant perturbation of the 
surface. A white noise is first generated using Phillips' spec- 
trum, then frequencies are divided in two sets. Low frequen- 
cies accounting for the global motion of waves are stored 
as a displacement map in a vertex shader. High frequencies, 
representing fine details, are stored in a normal map used 
for rendering. Since only a part of the spectrum is animated, 
a realistic ocean surface is obtained at interactive rates; an- 
other implementation by Chiu and Chang [CC06] combine 
this approach with an adaptive tessellation method [Joh04] 
(see previous section). 

An adaptive spectrum sampling is presented by Frechot 
[Fre06] in order to obtain a more realistic surface with capil- 
lary waves or riddles. The main idea of this approach is to re- 
sample the spectrum corresponding to high frequencies ar- 
eas, depending on the distance to the viewer. Robine and Fre- 
chot [RF06] also described a real-time additive sound syn- 
thesis method applied to the ocean surface. Sound synthesis 
methods are commonly used to efficiently sum a large num- 
ber of sinusoidal components called partials. Here, ocean 
waves are considered as partials by shifting from the fre- 
quency and time domains to the wave number and spatial 
ones. 

2.3. Hybrid approaches 

Hybrid approaches are a combination of spatial and spec- 
tral definitions, generating a geometric representation of the 




Figure 5: Ocean surface obtained in [TDGOO] 



surface while describing the components of wave trains in a 
realistic manner. 

The method proposed by Thon et al [TDGOO, TG02] 
shows a combination of a linear superposition of trochoids, 
which characteristics are synthesized using Pierson and 
Moskowitz's spectrum. In order to get fine scale details, a 
3D turbulence function is added to the surface; this method 
provides very realistic results, as shown on Figure 5. 

Lee et al. [LBR07] combine a superposition of sinusoids 
and the TMA spectrum (Texel, Marson and Arsole, based on 
Jonswap's) [BGRV85], creating a better statistical distribu- 
tion of waves. It also offers better control for the end-user 
who can change parameters such as the bottom depth, thus 
enabling him to generate different types of water volumes 
(lakes, rivers, etc). 

Xin et al [XFS06] extend the method of [TG02] to simu- 
late the effects of wind variations over the surface by modu- 
lating each spectrum frequency accordingly. 

Hinsingert et al [HNC02] apply a LOD scheme to the 
method of [TDGOO] by reducing both the sampling resolu- 
tion of the surface mesh and the number of trochoidal com- 
ponents depending on the distance to the viewer. This ap- 
proach reduces the amount of computations for the ocean 
simulation and tends towards real-time rates (approximately 
20 to 30 frames per sec. on a GeForce 2 GPU), by using only 
a subset of the components. 

A real-time simulation method for ocean scenes was pro- 
posed by Lachman [Lac07] with a framework represent- 
ing the surface either in the spatial domain by combining 
a spectrum and a sum of sinusoids or by using the approach 
of [TesOl]. This method is highly controllable, allowing the 
user to modify parameters such as the bottom depth or the 
agitation of the sea (defined by the Beaufort Scale). A LOD 
representation of the surface is also implemented to reduce 
computation costs. The main idea consists in defining sev- 
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eral resolution levels from the projection of the pyramid of 
vision on the surface, and meshing the surface accordingly. 

2.4. Discussion 

The main advantage of spatial domain approaches is their 
ability to produce a simple and fast simulation of the ocean 
surface, but they require a large number of periodical func- 
tions for visually plausible results. Several optimizations 
have been proposed such as GPU evaluation or adaptive 
schemes to reduce this number according to the distance to 
the viewer. However, using sine or cosine functions induces 
a too rounded shape of waves and the surface appears too 
smooth. 

Spectral domain methods address this problem by using 
oceanographic data directly and allow to disturb the surface 
according to physical parameters such as wind speed, which 
brings more realistic results. Unfortunately these methods 
lack from global control of ocean waves. 
In this context, hybrid approaches represent a good compro- 
mise between the visually plausible results obtained by spec- 
tral approaches and the global control offered by spatial do- 
main methods, and allow to obtain fast and effective simula- 
tions of the ocean surface's dynamic behaviour. A complete 
summary of the methods presented in this section is shown 
on Table 1 . 

However, all these methods only consider deep water phe- 
nomena, in which the ocean surface is being subjected to 
small perturbations. Indeed, in shallow water parametric or 
spectral approaches cannot faithfully reproduce ocean dy- 
namics near coasts, e.g. breaking waves. To address this 
complexity, we have to focus on approaches that consider 
interactions and collisions between the ocean and the shore, 
which are presented in the next section. 

3. Ocean dynamics simulation in shallow water 

Navier-Stokes equations (NSE) are able to capture the dy- 
namic complexity of a fluid flow. Incompressible flow is usu- 
ally considered for large volumes of fluid like oceans: 

V(/=0 (11) 

f + i/V(; + ^-^,^-g = o (12) 

with U = {u,v,w) the velocity of the fluid, /j its viscos- 
ity, p its pressure, p its density, and g representing gravity 
(0,9.81,0). Operator VU (resp. V^(7) is the gradient (resp. 
Laplacian) of U and is defined by VO = ^, ^) (resp. 

V^t) = 1^ + 1^ + 1^). The first equation guarantees mass 
conservation, i.e. the density of the fluid remains constant 
over time. The second one guarantees moment conservation: 
the acceleration ^ of the fluid is equal to the sum of the 
applied forces weighted by its density. Other formulations 
for fluid behavior include shallow water equations (or Saint- 
Venant equations), preferably used whenever the horizontal 
scale is greater than the vertical one. 



There are two main ways of discretizing these equations. 
On the one hand, Eulerian approaches use a 2D or a 3D 
grid where the fluid goes from cell to cell. On the other 
hand, Lagrangian approaches are based on particles carry- 
ing the fluid, thus representing small fluid volumes. Hybrid 
techniques were also developed to combine these two tech- 
niques, adding small-scale details to Eulerian simulations 
such as bubbles or spray. An excellent tutorial on physically- 
based simulation of fluids for Computer Graphics is pre- 
sented by Bridson et al [BMFG06]. Adabala and Manohar 
[AM02] proposed a survey of most physically-based fluid 
simulation techniques. Since numerous phenomena can be 
simulated by these methods (fire, smoke, lava, etc.) we will 
focus on liquids and more specifically on ocean waters. 



3.1. Eulerian approaches 

Kass and Miller [KM90] were the first to use a physically- 
based realistic approach to represent the ocean surface. This 
method consists in solving Saint- Venant equations in 2D in 
order to obtain a heightmap. This was later extended by 
Chen and Lobo [CL95] to solve NSE taking into account 
pressure effects, hence generating a more realistic surface 
modulation. 

A full 3D resolution of NSE was introduced in the Com- 
puter Graphics community by Foster and Metaxas [FM96, 
FM97], inspired from a physical approach proposed by Har- 
low and Welch [HW65] based on a uniform voxel decom- 
position of the 3D space. Velocity of the fluid is defined at 
the centers of each face of voxels; gradient and Laplacian 
are computed using finite differences between faces of ad- 
jacent voxels. Divergence at a given voxel is obtained from 
the updated velocities at its faces. In order to guarantee fluid 
incompressibility, the pressure field at a given voxel is modi- 
fied according to the divergence. Virtual particles are placed 
in the voxel grid to extract the fluid surface, i.e. to discover 
its location, however the resulting surface does not look visu- 
ally convincing. This problem was later addressed by Foster 
and Fedkiw [FFOl] who represent it as an implicit, level- set 
surface evolving according to the fluid's velocity. Neverthe- 
less, visual results show compressibility artifacts since a no- 
ticeable volume of fluid disappears during the animation. 

In order to solve this problem, Enright et al [EMF02] 
place particles above and below the surface, and advect 
them according to the velocity of the voxel they belong 
to. The surface is then computed by interpolating between 
the two implicit functions obtained from the locations of 
these particles. Since this approach called particles level- 
set method can be applied to all types of fluids, it became 
very popular and was widely used for many purposes, in- 
cluding the simulation of interactions between fluids and 
solids [CMT04, WMT05] or between different types of flu- 
ids [LSSF06]. The main drawback is high computational 
costs, which can be reduced by using an adaptive, octree- 
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Model 


Implementation 


Spatial 


Spectral 


Wave refraction 


GPU 


Level-Of-Details 


1980 


Schachter [Sch80] 


X 










1981 


Max [MaxSl] 


X 










1986 


Fournier et al. [FR86] 


X 










1987 


Mastin et al. [MWM87] 




X 








1987 


Ts'oetal. [TB87] 


X 




X 






1999 


Gonzato et al. [GS99] 


X 










2000 


Gonzato et al. [GSOO] 


X 




X 






2000 


Thon et al. [TDGOO] 


X 


X 








2001 


Premoze et al. [PAOO] 




X 








2001 


Tessendorf [TesOl] 




X 








2001 


Schneider et al. [SWOl] 


X 






X 




2002 


Gamito et al. [GM02] 


X 




X 






2002 


Hinsinger et al. [HNC02] 


X 


X 




X 


X 


2002 


Isidore et al. [IVB02] 


X 






X 




2002 


Thon et al. [TG02] 


X 


X 








2003 


Cieutat et al. [CGG03] 




X 


X 






2004 


Finch [Fin04] 


X 






X 


X 


2004 


Johanson [Joh04] 


X 






X 


X 


2005 


Kryachko [Kry05] 


X 






X 


X 


2005 


Mitchell [Mit05] 




X 




X 




2006 


Chiu et al. [CC06] 




X 




X 


X 


2006 


Frechot [Fre06] 




X 




X 


X 


2006 


Huet al. [HVT*06] 




X 




X 


X 


2006 


Lee et al. [LGL06a] 


X 






X 


X 


2006 


Lee et al. [LGL06b] 


X 






X 


X 


2006 


Robine et al. [RF06] 




X 




X 




2006 


Xin et al. [XFS06] 


X 


X 




X 


X 


2007 


Chen et al. [CLW07] 


X 






X 




2007 


Chou et al. [CF07] 


X 




X 


X 




2007 


Lachman [Lac07] 


X 


X 




X 


X 


2007 


Lee et al. [LBR07] 


X 


X 




X 




2007 


Salgado et al. [SC07] 


X 






X 





Table 1: Classification of the simulation methods for deep water presented in section 2 



like structure [LGF04] to refine the voxel grid only near the 
fluid-solid or fluid-fluid interface. 

This approach permits to capture the dynamic behavior of 
fluid surfaces in a realistic way and can be used for oceanic 
scenes. Enright et al [EMF02] simulate breaking waves by 
constraining the velocities in the grid using parametric equa- 
tions [R098]. Although realistic results are obtained, this 
approach can only handle limited types of waves (spilling 
or rolling waves). Mihalef et al [MMS04] simulate break- 
ing waves by constraining velocities inside the grid using 
the parametric approach of Thon et al [TDGOO] in 2D. The 
initial horizontal and vertical velocities of each element of 
fluid are given by: 

{u =Awe^- co&{kx) 
V = Awe*' sin(^t) 



where, as in [TDGOO], A is the amplitude of a wave, k = 
{kx,kz) its characteristic vector and w its pulsation. The user 
can obtain any type of shore-break by choosing from a li- 
brary of precomputed profiles. The full 3D simulation is then 
computed by extruding the desired profile along the parallel 
direction to the shore (see Figure 6). 



Thiirey et al [TMFSG07] simulate breaking waves in real- 
time, by solving Saint- Venant equations in 2D. A breaking 
wave is obtained by detecting regions with high curvatures 
on the heightmap representing the fluid, and adding a new 
mesh to account for the water falling from the crest. Each 
vertex of this mesh is projected in the 2D grid and advected 
using the corresponding cell. 
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Figure 6: Eulerian breaking waves simulation from 
[MMS04] 



3.2. Lagrangian approaches 

In Lagrangian approaches, the fluid is represented by a set of 
particles following physical laws. Particle systems were first 
introduced in the computer graphics community by Reeves 
[Ree83] to simulate natural phenomena such as water, fire, 
clouds or smoke. Miller and Pearce [MP89] simulate in- 
teractions between particles by connecting them with the 
aid of springs to represent attraction or repulsion forces be- 
tween neighboring particles, hence simulating viscous liquid 
flows. Terzopoulos et al [TPF89] later extended this method 
to simulate morphing from a solid to a liquid by modifying 
the stiffness of each spring. Tonnesen [Ton91] simulates the 
same kind of effects by solving heat transfer equations on 
the particles set. 

Although realistic results are obtained with these meth- 
ods, they are unable to simulate the inner dynamic com- 
plexity of fluids described by NSE. Stam and Fiume [SF95] 
introduced the Smoothed Particle Hydrodynamics (SPH) 
model in computer graphics, originally used in astrophysics 
to study the dynamics of celestial objects as well as their in- 
teractions [Luc77], and formalized by Monaghan [Mon92] 
for fluids. In that case, NSE are simplified and linearized. 
By considering that each particle has its own mass and that 
their sum gives the overall mass of the fluid. Equation 11 
can be omitted since the overall mass remains constant over 
time. In equation 12, the non-linear advection term, repre- 
senting the fluid's transport related to its velocity in the Eu- 
lerian case, can be omitted: the fluid is transported by the 
particles themselves. Finally only one linear equation is nec- 
essary to describe the fluid's velocity U : 

^ = i(-Vp+^v2(/ + pg) (13) 
at p 

Let / = — Vp + /jV^t/ + pG represent the total amount of 
forces applied on a particle. The term V p accounts for pres- 
sure forces, jjV^U for viscosity forces, and pG for exterior 



forces. The acceleration a of a particle is then given by: 

" = ? = -(^'""' + /"'■" + /"') (14) 
at p ' 

where fP''"^^ (resp. /"^^ and f"'^') represents pressure (resp. 
viscosity and exterior forces). The main idea in the SPH ap- 
proach consists in defining particles as potential implicit sur- 
faces (also known as blobs or metaballs) influencing their 
neighbors. A particle represents an estimation of the fluid in 
a given region rather than a single molecule. Force compu- 
tations rely on the interpolation of any scalar value 5, for a 
given particle depending on its neighbors: 

Si=Y,mj-!-W{r,R) (15) 
;=1 fj 

where Np is the total number of particles, mj (resp. Sj) is the 
mass (resp. the scalar value) of particle j, r is the distance be- 
tween particles i and j, R is the maximal interaction radius, 
and W is a symmetric interpolation function. The gradient 
VSi and Laplacian V^5; are given by: 

N g . 

V5,- = mj-!-\/W{r,R) (16) 
j=l Pi 



c 

V^^i = 52 mj-!-\7^W{nR) (17) 
j=i Pj 

Equation 15 can then be used to compute the density of a 
particle /; pressure and viscosity computation will rely on 
equations 16 and 17. 

The SPH model was extended by several authors [DC96, 
SAC*99] to represent elastic deformations, lava flows, etc. 
by simulating heat transfer and varying viscosity. Muller 
et al [MCG03] presented an extension to simulate low 
compressible liquids such as water. Clavet et al [CBP05] 
represent visco-elastic fluids by coupling particles with 
a mass-spring system; particles advection is governed by 
Navier-Stokes equations but also Newtonian spring dynam- 
ics, which slows the fluid down depending on springs stiff- 
ness. Another extension by Muller et al [MSKG05] rep- 
resent interactions between fluids with different phases by 
modifying the viscosity between neighboring particles. This 
method can also simulate bubbles generated by heating, 
using varying densities for fluid particles to represent the 
amount of gasification. 

However, two problems are inherent to the SPH approach. 
First of all, compressibility effects can be observed when the 
volume does not remain constant over time, yielding non- 
realistic visualizations. Besides, a huge amount of particles 
is needed to simulate large volumes of water with fine-scale 
details, leading to high computation and memory consump- 
tion. These problems were addressed by several authors. 

Premoze et al [PTB*03] introduced the Moving Parti- 
cle Semi-implicit method (MPS) in computer graphics from 
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Figure 7: Hybrid simulation from [LTKF08] 



physically-based models [K096, YK096]. This approach 
is an extension of SPH where the density of a particle is 
modified to guarantee constant pressure through the solv- 
ing of a Poisson equation. It was extended by Wang et al 
[WZC*06] to simulate 2D breaking waves, combined to ob- 
tain a 3D surface. Becker and Teschner [BT07] guarantee 
near-incompressibility by modifying pressure depending on 
density variations. Since pressure is directly related to attrac- 
tion or repulsion forces between particles, it is possible to 
ensure a near constant distance between neighboring parti- 
cles, hence a near constant volume. An alternative method 
was also presented more recently by Solenthaler and Pa- 
jarola [SP09]. 

Adaptive schemes were proposed to reduce the number of 
particles and limit computation and memory costs for sim- 
ulating large volumes of fluid. Desbrun and Cani [DC99] 
merge or split particles according to their density. The over- 
all mass of the system is kept constant by re-computing par- 
ticle masses as soon as they are subjected to one of these 
operations. In order to guarantee symmetric interactions be- 
tween particles with different masses, a shooting / gathering 
scheme ensures that attraction or repulsion forces are equal. 
Although this approach reduces computation costs, visual 
artifacts can be observed during animations because of in- 
stantaneous merging or splitting of particles near the sur- 
face. The adaptive scheme proposed by Hong et al [HHK08] 
consists in varying the size and number of particles accord- 
ing to their position relative to a fixed set of layers de- 
fined by the user in a preprocessing step. Finally, Adams et 
al [APKG07] obtain adaptively sized particles by defining 
merging or splitting processes, depending on the distance to 
the viewer or to the surface. This ensures that visual artifacts 
are avoided near the surface since only particles deep inside 
the fluid can be modified. This approach was implemented 
on GPU by Yan et al [YWH*09] to simulate breaking waves 
in real-time. 



3.3. Hybrid approaches 

The main idea in hybrid approaches is to simulate the main 
body of fluid using an Eulerian method, and fine-scale de- 
tails such as foam, spray or bubbles with a Lagrangian 
method. By adding small details, very realistic results are 
obtained, and interest for this approach is growing more and 
more with the development of efficient simulations. O'Brien 
and Hodgins [OH95] couple an Eulerian method with a par- 
ticle system to represent sprays. Navier-Stokes equations are 
solved in 2D, and sprays are generated depending on the 
magnitude of the vertical velocity at each cell; these parti- 
cles are then advected using the velocity in the grid. The 
Eulerian simulation of Takahashi et al [TFK*03] is coupled 
with particles generated in high curvature regions, then ad- 
vected independently. Greenwood and House [GH04] use a 
level-set approach [EMF02] coupled with particles gener- 
ated in high curvature regions, advected using the velocity 
of their corresponding cell. The friction of the fluid is also 
taken into account to generate bubbles. This method was later 
extended by Zheng et al [ZYP06] by deforming each bubble 
based on surface tension, thus avoiding non-realistic spheri- 
cal bubbles. 

Kim et al [KCC*06] extend an adaptive level-set ap- 
proach [LGF04] to generate particles according to the 
temporal variation of the volume of fluid in each cell. 
The main idea consists in using particles advected below 
the surface, transformed into water or air at the rendering 
stage. Since the number of transformed particles depends 
on the volume of water lost in a cell, this characterizes 
the turbulence in that cell and hence generates particles in 
highly perturbed regions. 

Yuksel et. al [YHK07] introduced the concept of wave 
particles to simulate waves and their interactions with solid 
objects. Wave particles are defined by a disturbance func- 
tion characterized by an amplitude, a propagation angle and 
a maximal interaction distance. The fluid itself is represented 
by a height-field whose values are computed according to 
neighbouring wave particles. When waves collide with solid 
objects, propagation angles are re-computed according to the 
curvature at collision points and new particles are generated 
to obtain waves refraction. 

Thiirey et al [TRS06] proposed to couple 2D and 3D sim- 
ulations using Saint- Venant equations solved with a Lattice- 
Boltzmann method [FDH*87]. The 2D simulation accounts 
for global motion of the surface, but fine-scale details are 
represented in 3D. These two layers interact with each other 
to generate foam and sprays, by creating particles accord- 
ing to the velocity in each cell and surface tension at the 
air-water interface. This work was later extended with SPH 
to represent interactions between particles in [TSS*07]. The 
same approach is followed by Losasso et al [LTKF08] to 
combine an Eulerian method with SPH particles whose num- 
ber is computed from the divergence of the fluid's veloc- 
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ity, thus simulating bubbles or sprays generated by breaking 
waves (see Figure 7). 

3.4. Discussion 

The main advantage of Eulerian approaches is their ability 
to simulate a large scope of different phenomena, explaining 
why they received much attention from the computer graph- 
ics community. Nevertheless, for fine-scale details such as 
spray or bubbles, NSE need to be finely discretized which in 
turn demands high computations and memory consumption. 
Another drawback is the use of fixed grids, forbidding the 
fluid to flow outside the grid. 

On the other hand, Lagrangian approaches can be used 
to simulate a wide range of phenomena. Their main advan- 
tage is their ability to represent fine-scale details and to flow 
anywhere in a virtual environment, but a large number of 
particles is usually needed to obtain realistic results. This 
problem can be alleviated using adaptive split-and-merge 
schemes to reduce computation costs. However it is worth 
noticing that most of the computation time for one particle is 
spent in testing neighboring particles or other objects for col- 
lision. Therefore a broad-phase collision detection is usually 
implemented by storing particles in a virtual grid, meaning 
that Eulerian or Lagrangian approaches share common prob- 
lems such as defining an appropriate size for the grid's cells. 
This is also illustrated by hybrid methods which produce re- 
alistic results by adding details to Eulerian approaches using 
particle systems. The literature presented in section 3 is sum- 
marized on Table 2. 

4. Realistic ocean surface rendering and ligliting 

In the previous sections we presented different approaches to 
obtain a plausible shape of the ocean surface, in deep water 
or near the shoreline. This surface can then be rendered using 
water shaders implementing light reflection and refraction 
described by Fresnel equations. Nevertheless, the ocean's vi- 
sual aspect is characterized by numerous interactions. Their 
impact is visible through several phenomena such as foam, 
sprays and light interactions, addressed by several methods 
in computer graphics. Some of them were already mentioned 
earlier, when the fluid simulation algorithm includes these 
phenomena directly {e.g. bubbles). However in most cases 
rendering is performed in a post-processing step. We focus 
on the following on different methods designed to enhance 
the realism of ocean scenes. 

4.1. Foam and spray 

Foam and spray rendering methods can be divided into two 
types: on the one hand, methods based on the computation of 
foam amount at one point on the ocean according to empiri- 
cal models, and on the other hand, methods based on particle 
systems and rendered as point clouds. 




Figure 8: Foam rendering from [JGOl ] 



4.1.1. Empirical models 

Jensen and Golias [JGOl] compute the amount of foam at 
a given point according to its height: if the height differ- 
ence with its neighbors is greater than a predefined thresh- 
old, foam is generated and rendered using a semi-transparent 
texture (see Figure 8). This method was extended by Jeschke 
et al [JBS03] to take into account the amplitude of ocean 
waves. The main drawback with this approach is encoun- 
tered when trying to animate foam, since its motion does 
not follow the waves. Another problem arises from the em- 
pirical models that do not consider atmospheric conditions 
which are, in reality, the main cause of the apparition of 
foam. Oceanographic observations have shown for example 
that foam only appears when wind velocity exceeds 13km/h 
[Mun47]. 

Premoze and Ashikhmin [PAOO] compute the amount of 
foam based on wind velocity and the temperature difference 
between water and air. The user can control atmospheric 
conditions to induce more or less foam. However, results are 
limited by several constraints on the surface, and the anima- 
tion is not visually convincing. Darles et al [DCG07] extend 
this method with an animated texture, combined with other 
phenomena such as light attenuation due to spray particles 
in the air. 

4.1.2. Particle systems 

Peachey [Pea86] first proposed to use particle systems in or- 
der to represent spray generated by breaking waves. Wang et 
al [WZC*06] extend this approach by describing breaking 
waves using a Lagrangian approach, and generating spray 
subsystems according to the velocity of water particles (see 
Figure 9). Holmberg and Wunsche [HW04] generate foam 
particles by considering the amplitude of the waves, and ren- 
der them as cloud points. The same idea is used by other 
authors [JGOl, JBS03,CC06] to generate spray particles ac- 
cording to the temporal variation of height at a given point. 
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Model 


Resolution 


Implementation 


Eulerian 


Lagrangian 


2D-based 


Full 3D 


Incompressibility 


LOD 


1989 


Miller et al. [MP89] 




X 




X 






1989 


Terzopoulos et al. [TPF89] 




X 




X 






1990 


Kass et al. [KM90] 


X 




X 








1991 


Tonnesen [Ton91] 




X 




X 






1995 


Stam et al. [SF95] 




X 




X 






1995 


Chen et al. [CL95] 


X 




X 








1995 


O'Brien et al. [OH95] 


X 


X 


X 








1996 


Desbrun et al. [DC96] 




X 




X 






1996 


Foster et al. [FM96] 


X 






X 






1997 


Foster et al. [FM97] 


X 






X 






1999 


Desbrun et al. [DC99] 




X 




X 




X 


1999 


Storaetal. [SAC*99] 




X 




X 






2001 


Foster et al. [FFOl] 


X 






X 






2002 


Enright et al. [EMF02] 


X 






X 


X 




2003 


Premozeet al. [PTB*03] 




X 




X 


X 




2003 


MuUer et al. [MCG03] 




X 




X 






2003 


Takahashi et al. [TFK*03] 


X 


X 




X 


X 




2004 


Carlson et al. [CMT04] 


X 






X 


X 




2004 


Greenwood et al. [GH04] 


X 


X 




X 






2004 


Losasso et al. [LGF04] 


X 






X 


X 


X 


2004 


Mihalef et al. [MMS04] 


X 




X 




X 




2005 


Clavet et al. [CBP05] 




X 




X 






2005 


Muller et al. [MSKG05] 




X 




X 






2005 


Wang et al. [WMT05] 


X 






X 


X 




2006 


Kim et al. [KCC*06] 


X 


X 




X 




X 


2006 


Losasso et al. [LSSF06] 


X 






X 


X 




2006 


Thurey et al. [TRS06] 


X 


X 


X 


X 


X 




2006 


Wang et al. [WZC*06] 




X 


X 


X 


X 




2006 


Zheng et al. [ZYP06] 


X 


X 




X 


X 




2007 


Adams et al. [APKG07] 




X 




X 




X 


2007 


Becker et al. [BT07] 




X 




X 


X 




2007 


Thurey et al. [TMFSG07] 


X 




X 


X 






2007 


Thurey et al. [TSS*07] 


X 


X 


X 




X 




2007 


Yuksel et al. [YHK07] 


X 


X 


X 


X 






2008 


Hong et al. [HHK08] 




X 




X 




X 


2008 


Losasso et al. [LTKF08] 


X 


X 




X 


X 




2009 


Solenthaler et al. [SP09] 




X 




X 


X 




2009 


Yanetal. [YWH*09] 




X 




X 




X 



Table 2: Classification of the simulation methods for shallow water presented in section 3 
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Figure 9: Spray rendering from [WZC*06] 



Although these methods provide realistic results, they also 
demand high computations since a large number of parti- 
cles is needed; this proves difficult to handle when rendering 
large oceanic scenes. 

4.2. Light-water interactions 

Light-water interactions play an important role in the way 
we perceive ocean's color, however they must be approxi- 
mated especially for real-time implementations. We use their 
order of approximation to classify existing methods, from 
strong simplifications of optics laws to highly accurate com- 
putations generating rippling caustics formed when light 
shines through ocean waves. 

4.2.1. First order approximation 

A first approach consists in rendering light-water interac- 
tion by applying geometrical optics laws, using Fresnel co- 
efficients or Schlick's approximation [Sch94]. The Beer- 
Lambert law reduces the intensity / of a light ray by an ex- 
ponential factor, depending on the covered distance d (i.e. 
the distance traveled by the ray from the light source to the 
intersection point) and an attenuation coefficient a: 

-ad 



l = l{0)e~ 



(18) 



where /(O) is the initial intensity of the ray. In the case of 
reflexion, intensity can be reduced according to the optical 
depth due to particles in the air. For refraction, the reduction 
comes from the attenuation coefficient of water, as imple- 
mented by Tessendorf [TesOl]. 

Gonzato and Le Saec [GSOO] use Ivanoff's table defin- 
ing the color of an ocean depending on its position on Earth 
[JN74]. This color takes into account reflected and refracted 
light, and in presence of submarine objects, a fast light ab- 
sorption method is applied by using an ocean water absorp- 
tion spectrum. 




Figure 10: Creation of illumination volumes from [IDN02] 



In order to reduce computation costs and obtain realistic 
results in real-time, simpler methods and GPU implementa- 
tions are often preferred. An environment mapping method 
can be enough to approximate reflexion and refraction ef- 
fects. In both cases, the viewer is assumed to be far from 
the considered point, thus incident rays are considered par- 
allels and reflected rays only depend on the surface's normal 
vector. This technique is commonly used in computer graph- 
ics, especially for ocean surface rendering [SWOl, JGOl]. 
Baboud and Decoret [BD06] represent refraction by an envi- 
ronment map of the bottom of the sea. An alternative method 
was proposed by different authors [Bel03,HVT*06] by map- 
ping a sky texture onto the surface and considering this sur- 
face as perfectly specular. 

Recently Bruneton et al. [BNHIO] proposed a GPU im- 
plementation of ocean lighting based on a hierarchical rep- 
resentation mixing geometry, normals and BRDF. The re- 
flected light coming from the sun and the sky and the light re- 
fracted by the ocean are simulated by combining the BRDF 
model of Ross et al. [VRP05] and wave slope variance (i.e. 
normals) depending on the required level of detail. 

4.2.2. Multiple order approximation 

When a light ray enters the ocean surface, a fraction is dif- 
fused, absorbed and reflected on organic particles in suspen- 
sion near the surface. In order to represent the complexity 
of this phenomenon and compute the amount of light at any 
point inside the water volume, a Radiative Transfer equation 
(RTF) must be solved. 

Arvo [Arv86] propose to emit rays from light sources and 
to accumulate them in illumination maps to represent caus- 
tics. Iwasaki et al [IDNOl, IDN03] use beam-tracing meth- 
ods [HH84, Wat90] to solve RTF by discretizing the water 
volume into illumination volumes formed by the intersec- 
tion of refracted rays and the bottom of the sea, sliced into 
sub-volumes (see Figure 10). The amount of light bounc- 
ing back to the surface is computed by accumulating light in 
all traversed sub-volumes. By discretizing the water volume 
into horizontal, planar sheets, bouncing light can be com- 
puted by accumulating light at intersection points between 
rays and each sheet to obtain caustics. The photon mapping 
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Figure 11: Ocean surface rendering for a clear water from 
[PAOO] 



method proposed by Jensen [JenOl] can also simulate caus- 
tics or other global illumination effects. 

Nishita etal [NSTN93] proposed a rendering algorithm to 
visualize oceans viewed from space by taking into account 
optical properties due to atmospheric conditions. The color 
of the ocean is computed using a spectral model based on the 
RTE which takes several physicals parameters into account 
such as ocean depth, incident angles of the sun and viewing 
direction. 

Sub-surface scattering was also addressed by several au- 
thors [PAOO, CS04] by taking into account several bio- 
optical parameters such as chlorophyll concentration and 
turbidity [Jer76]. In that case, coefficients for attenuation 
a{X) and diffusion b{X) are obtained from physical laws 
[Mor91,GM83] according to the phytoplankton concentra- 
tion Cp inside water: 

a{X) = (fl„,(?i)+0.06C°-''^)(l+0.02e-°'""*(^-^'^"') 



where Ow (A.) is the attenuation coefficient of pure water for 
a given wavelength X. Thus turbidity is directly related to 
the concentration of phytoplankton. This approach is able to 
represent realistic ocean colors according to their biological 
properties (see Figure 11). 

Gutierrez et al [GSAM08] simulate a wider range of phe- 
nomena such as inelastic diffusion of light rays, fluorescence 
effects due to phytoplankton particles and Raman diffusion, 
observed when an incident ray's wavelength is modified be- 
cause of water properties. This model is able to character- 
ize the visual impact of different kinds of sediments on the 
resulting image, and thus represent different ocean waters 
depending on their bio-chemical composition. 

Finally another type of methods consider subsurface scat- 



tering inside water as participating media. Interested readers 
should refer to [PCPS97] for an overview of these methods. 

4.3. Discussion 

Numerous methods are available in the literature to take dif- 
ferent optical phenomena into account and increase the real- 
ism of oceanic scenes. In the case of foam and sprays, em- 
pirical methods can effectively simulate these phenomena 
according to the state of disturbance of the surface. Com- 
bined with particle systems, these approaches allow to ob- 
tain realistic details, however they require a large number of 
particles which induces an important memory cost and com- 
putation time at rendering stage. Moreover, empirical meth- 
ods are usually efficient for deep water scenes only; to our 
knowledge there is no work dealing with realistic foam cre- 
ated by breaking waves. 

For light-water interactions, different methods were pro- 
posed to represent the complexity of these non-linear phe- 
nomena, yielding more and more realistic results. Moreover, 
recent approaches even consider physical and biology effects 
to simulate a wide variety of virtual environments. All these 
works are summarized on Table 3. 

5. Conclusion 

In this survey we have presented available methods in com- 
puter graphics for modeling and rendering oceanic scenes. 
The great amount of different works is inherent to the ex- 
treme diversity of the phenomena involved. 

The first two sections of our survey focused on computer 
graphics models able to simulate the dynamical behavior of 
the ocean. We have seen that those methods are divided into 
two categories: on the one hand, methods usually dedicated 
to deep-water simulation, on the other hand fluid-based ap- 
proaches trying to represent breaking waves near the shore. 
In the last section, we have seen different methods able to 
represent several phenomena involved in ocean rendering, 
namely foam, sprays and light-water interactions, that make 
for visual realism. 

The next decade could see new methods emerging in an 
attempt to bridge the gap between deep-sea simulations and 
breaking waves representations; it would be necessary to put 
together different types of simulations on different scales. 
Scalable methods are also required for modeling and ren- 
dering stages in order to obtain real-time rates, which could 
include dynamic sampling of the simulation domain (as pro- 
posed by Yu et al [YNBH09] for rivers) and adaptive render- 
ing models taking the visual impact of the considered phe- 
nomena into account. 
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Foam and Spray 


Light-water interactions 


Empirical 


Particles 


First order 


Multiple order 


1986 


Arvo [Arv86] 








X 


1986 


Peachey [Pea86] 




X 






1993 


Nishita et al. [NSTN93] 






X 




2000 


Gonzato et al. [GSOO] 






X 




2000 


Premoze et al. [PAOO] 


X 






X 


2001 


Iwasaki et al. [IDNOl] 








X 


2001 


Jensen [JenOl] 








X 


2001 


Jensen et al. [JGOl] 


X 


X 


X 




2001 


Schneider et al. [SWOl] 






X 




2001 


Tessendorf [TesOl] 






X 




2003 


Belyaev [Bel03] 






X 




2003 


Iwasaki et al. [IDN03] 








X 


2003 


Jeschke et al. [JBS03] 


X 


X 






2004 


Cerezo et al. [CS04] 








X 


2004 


Holmberg et al. [HW04] 




X 






2006 


Chiu et al. [CC06] 




X 






2006 


Baboud et al. [BD06] 






X 




2006 


Hu et al. [HVT*06] 






X 




2006 


Wang et al. [WZC*06] 




X 






2007 


Darles et al. [DCG07] 


X 




X 




2008 


GutieiTez et al. [GSAM08] 








X 


2010 


Bruneton et al. [BNHIO] 






X 





Table 3: Classification of the realistic rendering methods presented in section 4 
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